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ABSTRACT 

We describe a accurate and fast pixel-based statistical method to interpolate fields of arbitrary spin 
on the sphere. We call this method Fast and Lean Interpolation on the Sphere (Flints). The method 
predicts the optimal interpolated values based on the theory of isotropic Gaussian random fields and 
provides an accurate error estimate at no additional cost. We use this method to compute lensed 
Cosmic Microwave Background (CMB) maps precisely and quickly, achieving a relative precision of 
0.02% at a HealPix resolution of = 4,096, for a bandlimit of ^ max = 4,096 in the same time it 
takes to simulate the original, unlensed CMB map. The method is suitable for efficient, distributed 
memory parallelization. The power spectra of our lensed maps are accurate to better than 0.5% at 
I = 3,000 for the temperature, the E and B mode of the polarization. As expected theoretically, we 
demonstrate that, on realistic cases, this method is between two to three orders of magnitude more 
precise than other known interpolation methods for the same computational cost. 



1. INTRODUCTION 

Gravitational lensing of the Cosmic Microwave Back- 
ground (CMB) is a unique probe of the distribution of 
mass in the entire visible Universe. The first signa- 
tures of this CMB lensing signal have already been seen 
in cross-correlation with large scale structure templates 
([Smith et alJl200ft iHirata et al1l2008D . The effect of lens- 
ing on the CMB power spectra will provide powerful ad- 
ditional constraints on the physics of the dark sector 
(jAbazaiian fc Dodelsonl 12003 : iSmith et all l2006afl . In- 
cluding reconstruction of the lensing potential in cosmo- 
logical parameter analysis removes further degeneracies. 
Once the lensing deflection is mapped with high preci- 
sion it will lik ely allow detecting the absolute mass scale 
of ne utrinos (Kaplin ghat et all 120031 : iLesgourgues et all 
120061 ) and provide tight constraints on the presence of 
dark energy at redshifts z ~ 2. Lensing creates po- 
larization B modes by rot ating the stronger E modes 
(jZaldarriaga fc Sel jak 1997]). These lensing B modes are 
not affected by cosmic variance and can therefore signif- 
icantly improve reco nstructions of the lensing potential 
([Smith et al.ll2006bh . The modeling of lensing B modes 
is further motivated since they are a foreground for the 
search for the ff-mode signal of inflationary gr aviational 
waves (|Knox fc Sondl2002l : iKesden et al1l2002h . 

Planck will yield the first all-sky temperature and 
polarization maps of the CMB with sufficient signal to 
noise to detect th is signal. Ground based expe riments 
such as QUBIC ( Hamilt o^iirCharlass^[2010h. SPT- 
pol ([McMahon et al.ll2009D or QUIET (jNewburgh et al.l 
120051 ) contain detailed information about the fine-scale 
structure imprinted by lensing while still covering signif- 
icant portions of the sky. For the extraction of lensing 
and polarization science from these data sets, fast and 
precise methods to simulated lensed all-sky CMB maps 
are indispensable. 



In the Born approximation, which specifies that the 
lensing effects may be modeled by the impact of a single 
lens plane on the CMB, the lensing corresponds to look- 
ing up at the value of the temperature or the polarization 
at displaced position on the sky. Thus, making lensed 
maps is effectively a resampling of the CMB on a dif- 
ferent set of positions. Different technical solutions were 
proposed. Currently three main techniques are known. 

The first technique consists in doing a brute force re- 
summation of the spherical harmonics at a new set of 
positions. Acceleration methods using symmetry prop- 
erties of the sampling positions can not be used. While 
prohibitively expensive for practical use, this method is 
exact and we therefore employ it on a subset of pixels as 
a precision benchmark for other methods. 

The second technique is polynomial interpolation on 
the sphere to compute the value of the temperature and 
polarization field at the displaced positions. The most 
simple of the interpolation technique is a bilinear in- 
terpolation of fields on the sphere (hereafter called the 
naive technique). Another algorithm, a Lagrange poly- 
nomial interpolation on the Equi-cylindr ical Projection, 
was first proposed by IHirata et al.l ([20041 ) and then later 
used by iDas fc Bodel (|2008l ) for a full-sky simulation of 
the CMB lensing on a light cone of a cosmological sim- 
ulation. A variant of this method, involving a bicu- 
bic interpolation scheme, is implemented in th e pub- 
licly available code LensPdQ d escribed in iLewisI (|2005[ ) 
and lHamimeche fc Lewlsl (|2008f ) (hereafter the ECP al- 
gorithm). 

A third t e chniq ue has been recently developed by 
iBasak et al.l (|2009l ). This algorithm consists in recast- 
ing the spherical harmonic coefficient of the unlensed 
field into the Fourier basis in the (#, (j)) variables. Then, 

1 http : //cosmologist . inf o/lenspix 
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we may use a Non-equispaced Fast Fourier Transform 
(NFFT) to compute the field at the displaced positions. 
This method achieves high precision and is significantly 
faster than the first technique though still too slow to 
allow for the production of a large set of lensed maps, 
which is required for the statistical analysis of observed 
CMB data. 

We propose a fourth technique that relies on the sta- 
tistical properties of the considered fields to be lensed. 
It is based on the idea of interpolating the original field 
but using the known spectral information to compute 
the correct weighing coefficients for the interpolation . 
This method is related to a Wiener filter ( Wiener|[l949[ ) 
but not limited to pure Gaussian random field. We call 
this method Fast and Lean Interpolation on the Sphere 
(Flints) In this work, we analyze the result obtained 
based on a HealPix pixelization. We note that our 
framework may be used on any pixelization of the sphere, 
including the Equi-cylindrical projection already used 
in previous works. The implementation that we pro- 
pose takes advantage of the geometrical properties of 
HealPix for a number of memory and computational 
optimizations. 

We note that our interpolation method is of general 
interest beyond its application to lensing. In fact, since 
it is based on Wiener filtering it is guaranteed to give the 
best possible mean squared error of any method for a field 
with the same power spectrum and for fixed interpolation 
stencil. 

The structure of this paper is the following. In Sec- 
tion [2j we describe the general interpolation method. 
Then, in Section [3j we discuss its performance in gen- 
erating lensed CMB maps. In Section |4j we conclude. 



2. INTERPOLATING CMB FLUCTUATIONS 

We propose a direct, simple though sufficiently precise 
method of interpolating complex fields on the sphere. 
This method is based on the most likely value an isotropic 
Gaussian random field takes at an arbitrary location on 
the sphere, given a set of sampled values (the interpola- 
tion stencil) and its power spectrum C\. In Section [2TT| 
we study the general method for the special case of spin- 
fields. Then, we generalize the obtained equations to 
spin-s fields in Section [2^21 We then detail in Section [231 
the algorithmic steps used to compute the identifiers of 
the pixels that are used in the interpolation. In Sec- 
tion ^. 4| we describe the memory and computational time 
optimizations that we used in Flints. We then list all 
the steps required to achieve the interpolation in Sec- 
tion 12.51 Finally, in Section 12. 6| we discuss the results of 
the raw performance tests of our interpolation method. 



2.1. Interpolating a scalar field 

We start by considering a Gaussian random scalar field 
on the sphere given by the function T(n), in the direction 
h. The general form of the joint probability of the value 

2 The reference implementation is written in C++/OpenMP and 
is available at http://www.iap.fr/users/lavaux/flints.php 



of T(n) in N + 1 directions is: 
P(T(n ),T(n 1 ),...,T(n N )) 



with W the inverse of the correlation matrix being 

W = V~\ (2) 



with Vij defined for two directions hi and hj as 



(3) 



with the correlation function £(cos0) for two directions 
separated by an angle 
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and the intrinsic variance of the temperature field 

<x 2 = C(cos(0) = 1). (5) 

Now, we can compute the conditional probability of the 
value of T(ho) given the N other values: 

V(T(h )\T(h 1 ),...,T(h N )) = 



^—exp I ^- (T(n ) —T) J , (6) 



with 



T=W f 



N 

o,o 
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and Wo,o the top-left most element of the W matrix. So, 
given the values in the directions {ni, . . . , n^v}, we define 
the interpolated value in the direction no to be equal to 
the most likely value T as defined above. The precision 
of the interpolation is given by the amount of allowed 
fluctuation 1 / ^Wo,o- The advantage of this procedure 
is that it is flexible in the number of points we take into 
account for the interpolation. Additionally, it remains 
purely local and the complexity scales as 0(N 2 x AT p i x ), 
with N the number of of points to compute the interpo- 
lation and iVpi x the number of pixels in the map. This 
locality allows us to take full advantage of the parallelism 
offered by current multi-core CPUs and by distributed 
computing environments. An illustration of our interpo- 
lation procedure is given Fig. [T] for the case of interpola- 
tion stencil with 9 elements ("neighbours"). 

We further increase the speed of the interpolation by 
precomputing the covariance matrix linked to the pix- 
elization, that is the inverse of the matrix Si j defined 
by: 



Si' 



-^{T(hi)T(hj)) 



(8) 



for hi the direction corresponding to the center of a pixel 
of the map to be interpolated. Strictly speaking, we 
would need to recompute these matrices if the angu- 
lar power spectrum of temperature fluctuations changes. 
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But as the weights are continuous functions of the angu- 
lar power spectra, two relatively similar spectra should 
give the same weights. This alleviates the need of recom- 
puting these weights for any single change of the angular 
power spectra, at the potential cost of a small loss of pre- 
cision. Si j is in practice a 9 x 9 matrix as we have nine 
neighbors for any direction of interpolation. Once we 
have this matrix it is fast to determine the Wp .j by doing 
block matrix computation ([Press et al.lll992f ). First, we 
write the shape of the V matrix: 



V 



1 B T 
B S 



(9) 



We let the block matrix shape of its inverse W have the 
following shape 



W 



B S 



(10) 



and write that the product should make identity: 
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which leads us to the following equalities: 
~ 2 -{l-B^S- l B), 
~ X B, 
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S=^S~ 1 BB ] S- 1 . 



(11) 



(12) 
(13) 
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The matrix S gives the correlation between the different 
direction of the pixelization, B gives the correlation be- 
tween the sought interpolated direction and the direction 
of the pixelization, <j is the standard deviation of the in- 
terpolator in the interpolated direction. The equations 
(fT2j) and ([T3]) are computed as needed for each interpo- 
lated direction. The equation (fT4|) is not used as we do 
not need this part of the matrix. 

Finally, we may express the value of the temperature 
in the interpolated direction n: 



f(n) 



N 



! (T(h)T(h 3 )} 



T(fn), 



(15) 







for which the Gaussian variance of the error of this esti- 
mator is 



N 



^ 5r/(T(n)r(^))(r(n)r(n j )) • (16) 



u i,j=l 

The above Eq. ([T5|) is the strict equivalent of Eq. (0 . We 
now have both an interpolated value and an estimate of 
the interpolation error. Even though we used a Gaussian 
field theory to derive these estimators, the interpolation 
is also optimal in a least-squared sense for non-Gaussian 
field. 




Figure 1. Illustration of our interpolation procedure - We present 
here an illustration of the procedure. To get the value of the field 
interpolated at the position illustrated by the red star, we use the 
value of this field sampled on the grid given by the squares. In this 
setting only the nine nearest neighbors, with temperature value Tj, 
contribute to the value Tq. 



The tim e complexity of a H ealPix spherical harmonic 
transform (jGorski et al.ll2QQ5[ ) is 0(N S i^ e £^ nax ), with ^ max 
the number of i modes used for the spherical harmonic 
transform. The difference in complexity between the 
HealPix spherical harmonic transform and our method 
means that if the number of neighbors N is sufficiently 
small then our method would be faster to compute a com- 
plete CMB sky in any direction than generating a very 
high resolution HealPix map. This improvement may 
be of order 0(£^ iax /(A/' 2 7V S ide,high)) for large £ max , where 
^side,high is the high resolution map needed for the naive 
technique. 

We note that, in the high resolution regime, the scal- 
ing of this algorithm is like the one of the ECP technique 
(jHirata et all 120041 : lHamimeche fc Lewisl 12008). On the 
other hand, there are several advantages to our proce- 
dure. First, the method provides an accurate error esti- 
mate for no additional computational cost. The weights 
take into account the part of the signal which is not well 
represented by the pixelization. The ECP interpolation 
implicitly assume that all the information is stored com- 
pletely in a neighborhood of the direction of interpola- 
tion, whereas our interpolation takes into account the 
possible lack of complete information. Second, Flints 
does not need another specific spherical harmonic trans- 
form to generate the lensed fluctuation on the sphere. It 
may use a previously existing or independently generated 
map. 

2.2. Interpolating spins fields 

We may now adapt our method for the case of spin-s 
fields on the sphere. The concepts are the same, except 
that we must now handle complex fields. We want to 
interpolate the complex field P(h). We assume that this 
field is isotropic and its correlation function is defined 
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by: 
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C Ptl ( s Y lt _ s ((3,a))e- is Xl7) 



= £(^)^- s ,-^)e-* s(7+ 
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= .C(«,/3,7), (19) 

with a,yS,7 the Euler angles defining the rotation to 
transform n' into n and s Yj >m the spin s spherical har- 
monic (see Appendix lAl for the definition). 

The function d!_ s _ s is the Wigner d function. An algo- 
rithm to compute this function is reca lled in Append ix [Bl 
and h as already been detailed in iTrapani fc Navazal 
([20061 ). An illustration that defines more properly these 
angles is given in Fig. [2j 

As in Eq. (|3]) , we define the Hermitian correlation ma- 
trix for two directions hi and hj, < z, j < N: 



«<(o,o,o) ' 



(20) 



with ctij, fyj and 7^ the Euler angles defining the ro- 
tation to transform hi in hj. Similarly as in Eq. (|9j), 
we define the Hermitian matrix S corresponding to the 
correlation of the values taken by the HealPix pixels, 
which corresponds here to the indices 1 < i, j < N. 

Using some geometry on the sphere it is possible 
to compute the above angles. The relation between 
h f (6 f ,(f) f ), h(9,(/)) and (a, /?, 7) is given by: 
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(21) 
(22) 

(23) 
(24) 
(25) 
(26) 



The choice of the sign of f3 is arbitrary, but the sign of 
the other angles depend on this original choice. As we 
need to compute the inverse of the above trigonometric 
identities we have decided to use the inverse of the tan- 
gent function to improve numerical stability, taking care 
of the signs to recover the correct a and 7 angles in the 
range [— 7r;7r]. 

The probability distribution of the complex field P(ho) 
given the other value of t his f ield in the direction f^, 
i = 1 ... N, is as in Section 12.11 



P(n )-P\ 2 ), (27) 



V(P(n )\P(n 1 ),...,P(n N )) = 

W ofi 



Wq,q 
2tt 



-exp 



with 



N 



P=W -^W^P(n t ). 



(28) 




Figure 2. Euler angle convention - We represent here the con- 
ventions use d fo r the orientation and the value of the Euler angles 
used in Eq. (19} . 



The Equations ([T2]), ([T3J) and (Q3} are stiu van d but this 
time for complex matrices and vectors. In that case 
is the conjugated transpose of matrix A. The resulting 
equation for the interpolated value P(h), taken to be P 
is 



N 



^Pi^P^fn)) 



sC(0,0,0) 



P(hi), (29) 



with S as defined above in this Section. The variance of 
the interpolated value is 



a 2 P (h) = s C(0,0,0)x 



N 



1 (P (n)P*(n,))(P(n J )P*(n)) 
s((0, 0,0)4 



(30) 



In the following, we focus on the actual implementation 
of the interpolation algorithm which uses the above two 
equations to compute the interpolated value in any di- 
rection. 

2.3. Identifying neighbors 

To be fast, the interpolation procedure must rely on 
a limited number of pixels, and more particularly the 
pixels just in the immediate vicinity of the direction of 
interpolation, using those pixels which carry most of the 
information about the interpolated point. 

We discuss here a particular implementation to identify 
neighbors that relies on the HealPix framework. We 
use a method that relies on the use of the neighbors () 
function, which returns the immediate nine neighbors, 
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original A^^e target N s ^ e t Neighbors Precision Computational time Memory HealPix time 



(serial, minutes) (Megabytes) (serial, minutes) 
1,024 2,048 4,096 9 2 x lO" 2 0.8 440 2.7 
1,024 2,048 4,096 36 1 x 10~ 2 5.1 5,500 2.7 
2,048 4,096 4,096 9 3 x lO" 3 3.3 1,700 5.5 
2,048 4,096 4,096 36 3 x 10~ 4 19 22,000 5.5 
4,096 8,192 4,096 9 4 x 10~ 4 14 6,800 12 



Table 1 

Performance for maps supersampling 

We give here the computing performance of our method to interpolate a scalar field containing power up to a bandlimit of t = 4, 096. The 
computational time corresponds to the single processor wall time taken by the algorithm to interpolate the map at the given resolution 
from a map at iV side to a map at 2 x N s [^ e . The HealPix computational time is estimated on the spherical harmonic synthesis transform 
at the target N s ^ e resolution. The memory consumption gives the size of the pixelization cache in Random Access Memory. The precision 
corresponds to the square root of the average variance of the error predicted by FLINTS, this square root is divided by the standard 
deviation of the maps. The quoted serial times have been measured on an Intel Xeon E5410 2.33 GHz based computer. The standard 
deviation of the maps is 39 /iK/K. 



sorted geometrically, see Fig. [TJ The time complexity 
of neighbors () is constant both in the NESTED or RING 
mode of HealPix. 

A fast way of extending this algorithm to a higher num- 
ber of neighbors consists in using the NESTED mode of 
HealPix. Thi s exten sion has already been described in 
IWandelt et al.l ([199 8). We only consider here the neigh- 
bors symmetric according to the direction in which to 
interpolate in the NESTED tree sense. The neighbors 
at a level y > 1 in a map at resolution can be 

found using a pixelization at N s ^ e /y. Their number is 
then 9 x 4 2y_1 . We derive their identifiers in NESTED 
mode by computing the first neighbors using the func- 
tion neighbors () of HealPix at resolution N s ^ e /y. We 
shift the bits of these identifiers by 2(y — 1) and fill up the 
lower bits with all the possible combinations. This pro- 
cedure yields all the 9 x 4 y_1 pixel neighbors of a given 
direction. 

This procedure is more attractive given our computa- 
tional constraints than using the alternative procedure 
query_disc() : it allows a stable number of pixels per 
neighbors, their geometrical distribution according to the 
central direction is stable and insensitive to numerical 
rounding errors, which makes it easier to tabulate S -1 
and it is relatively fast to compute the list. For these 
reasons, we only use this procedure in the rest of this 
work. 

2.4. Tabulating S' 1 

To reduce the time needed to compute the interpola- 
tion, we tabulate the matrices S -1 for all possible central 
pixels (corresponding to T 5 in the Figure [1]). This means 
we would need to precompute and store a matrix for 
each pixel of the map to interpolate. We made use of the 
symmetries of the HealPix pixelization to reduce this 
cost. As there is a rotational invariance in the equatorial 
part of the pixelization, it suffices to store one pixel per 
ring. For the polar regions, it would be sufficient to store 
the weights for only one quarter of one of the HealPix 
base tiles, since there is an eight-fold rotational invari- 
ance and an additional north-south symmetry. Using 
all these symmetries would reduce the amount of mem- 
ory required to store the precomputed weights by nearly 
1/48. 

For our implementation we use a subset of these sym- 
metries since there is a trade-off between reduction of 
memory use and the implementation complexity for the 



required pixel permutations. In addition to the equa- 
torial symmery, we use the north-south symmetry for 
the case of the spin-0 field only. Memory use scales as 
^(^ngb x ^s 2 ide)> w ^h N ng \) the number of neighbors cho- 
sen for computing the interpolation. 

To compute S -1 itself from 5, we use the Cholesky 
decomposition of S = t LL and do a dire ct inversion 
of L. We used the algorithm described in iPress et all 
(|1992f ). As the matrix S becomes so ill-conditioned that 
the Cholesky decomposition fail for high resolution map, 
we introduce an additional term n to help at computing 
this decomposition when this term is required. So, in 
practice, we decompose S + nX instead of S. By con- 
struction, the value of n is much smaller than one. If n 
is too big, the interpolation error is larger than it should. 
If n is too small, the decomposition fails. We use for n a 
value given by 

^ = ^min, negative H~ ^precision (^1) 

where A m in 5ne gative is the lowest negative eigenvalue of S 
at the used precision, e prec i s i n is a quantity dependent 
on the machine precision used for doing the actual com- 
putation on S. For double-precision, we take 1.49 10 -8 , 
which is the square root of the smallest deviation from 
1.0 detectable in a double precision representation. 

2.5. The interpolation algorithm 

Finally, the interpolation algorithm consists in achiev- 
ing the following steps: 

1. We set the resolution of the map to interpolate at 
A^ s ide and the band width to 4nax- 

2. We start by tabulating according to f3 G [0,7r] the 
C(0,/?,0) function, defined in Eq. (fT9j) . 

3. For the HealPix pixelization at N s ^ e ^ we compute 
S using Eq . ([20]) , invert it and tabulate it according 
to Section 12 A\ This is the end of the preparation 
phase. 

4. For any direction 0), we compute the identifiers 
of the pixel neighbors according to Section 12.31 

5. We sum up the value of the field sampled at those 
pixels and weighed them according to Eq. ([29]) . 

6. If required by the user, we also compute the error 
in the interpolation using Eq. ([30]) . 
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Interpolation absolute error 




I 



0.08 



muK/K 



Figure 3. Spatial distribution of predicted interpolation error for the upsampling — We show the distribution of the error predicted by the 
algorithm for each interpolated pixel of the sky. We predict a map at iV s ide = 4, 096 from a map at N s i& e /2 = 2,048, with Cax = 4, 096. 
The structure of the error distribution reflects the features of the HealPix pixelization. 



In all the subsequent tests, we have used the above 
scheme to generate the interpolated values. 

2.6. Supersampling interpolation performance 

We test the interpolator by supersampling a CMB map 
to twice the original resolution. Please refer to the results 
in Table [T] and in Figure [3l In all these tests, we start 
from a map at a resolution of iV^e and predict its values 
at twice the resolution, 2 x 7V S id e - Note that all these 
test are done for fixed bandlimit i = 4, 096 while the 
resolution of the grid is varied. For the lowest starting 
resolution of = 1,024 the field is somewhat un- 

dersampled, since modes beyond t = 2iV S ide begin to be 
noticeably aliased on the HealPix grid. Owing to the 
hierarchical property of the HealPix grid, the higher 
resolution pixels tile the lower resolution pixels. The dis- 
tance from the members of the interpolation stencil is 
therefore one fourth or three fourth of the size of a pixel. 
Based on numerical experiments we estimate the worst 
case error for directions farthest away from the mem- 
bers of the interpolation stencil to be no more than 50% 
higher than shown in Figure [3l 

In Figure [3j the distribution of the interpolation error 
reflects the features of the HealPix pixelization. The er- 
ror is uniformly small in the equatorial regime and shows 
the symmetries of the pixelization in the northern and 
southern caps. 

We tabulate the supersampling precision and the time 
consumption in Table HJ For comparison we list the time 
required to compute a single HealPix transform from 
{ai m } space to pixel space. The two operations take 
roughly the same time. Note that the spherical harmonic 
transform has only been done one way, whereas the inter- 
polation starts from pixels and yields pixels. As expected 
the error decreases with increased resolution and number 
of neighbors. It is interesting to see that we may get a 
decrease of one magnitude in the error by changing the 
number of neighbors at the resolution iV S ide = 2,048. Do- 
ing the same exercise with N s ^ e = 4,096 yields a large 



number of degeneracies in the S matrices. That shows 
that we reach the level where the problem of interpola- 
tion is dominated by errors in the floaing point represen- 
tation, as expected since a HealPix map at resolution 
TVgide is able to encode wavenumbers up to t ~ 27V S id e - 
So, as our £ max is 4,096 here, a map at 7V S ide = 2,048 can 
be interpolated at high precision with a sufficient number 
of neighbors. 

3. APPLICATION TO CMB LENSING 

In this section, we focus on the use of our interpola- 
tion procedure for producing lensed maps of the CMB. 
In Section 13. 1\ we recall the basic lensing equation in the 
Born approximation and the notation. In Section 13. 2\ 
we compare the predictive performance of FLINTS to 
LensPix and to the naive technique in HealPix. Fi- 
nally, in Section [331 we compare the power spectra of the 
lensed temperature and polarization as obtained through 
FLINTS and using CAMB. 

3.1. Theory of lensing 

Gravitational lensing acts like a remapping of the CMB 
photons on the sky. Thus the temperature signal on the 
sky is given by 

T bserved(n) = T C M B (n + d(n)), (32) 

where T b se rved is the observed temperature of the CMB 
in the direction n, Tcmb is the unlensed primary signal, 
and d(n) is the deflection field. This relation is also true 
for the polarization field defined by 



P(n) = Q(n)+iU(n) 



(33) 



with Q and U the Stokes parameters, with P being a 
spin-2 field. 

The deflection field defines in what direction and by 
what angle the photons were deflected from their original 
position. In the local basis defined by the direction n, 
we may define the angle a(n) as 



d(n) oc cos(a)ii6>(n) + sin(a)u0(n) 



(34) 
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time 
(minutes) 


Interp. 
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(minutes) 
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Precision 
time Z/2 
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1,024 


8 


1 


8 


2% 


11% 


32 


5.0 


9 


5% 


24% 


2,048 


9 


6 


20 


0.3% 


2.65% 


-160-250* 


12 


20 


9% 


125% 


4,096 


16 


22 


51 


0.04% 


0.3% 


-640-1,000* 


34 


52 


7.5% 


127% 



Table 2 

Lensing performance 

NB: Numbers with * are extrapolated from table 3 of Basak et al. (20Q3). "Interp." stands for "Interpolation". "Init." stands for 
"Initialization" . We measure the time for producing a lensed map, temperature and polarization, from a random realization of the CMB 
fluctuations and of the lensing potential. The total time is the sum of the interpolation time, the time to make an unlensed healpix map, 
the time to compute the deflection map but not the initialization time given in the first column. The L2 (Loo) precision corresponds to 
the standard deviation (maximum absolute value) of the error distribution divided by the standard deviation of the simulated CMB map, 
which is 39/iK/K. For all cases we used l max = 4,096 and nine neighbors for FLINTS, and a interpolation factor of one for the ECP 
method in LensPix. The maps at N s ^ e = 1,024 were therefore undersampled. 
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Figure 4. Interpolated CMB vs True value - Distribution of the 
error in the value given by both naive interpolation on a HealPix 
mesh and our method for different HealPix resolution. For each 
line, we represented the distribution of the relative difference be- 
tween the actual interpolated value and the true value. The 
green, black and blue lines correspond to HealPix interpolation 
at 7V side = 2,048, 7V side = 4,096 and 7V side = 8,192 respectively. 
The red line corresponds to our method at N s ^ e = 2,048. 



with (n,u<9,U0) the local spherical orthonormal basis, 
with u<9 = dn/d0 and = dn/d<p. In this basis, the 
lensed direction n' may be written as 

n' = n + d(n) = cos(|d|)n + sin(|d|) cos(a)ii0 

+ sin(|d|) sin(a)u0. (35) 

We prefer this form of the lensed direction instead of us- 
ing the angles for numerical stability, at the cost of com- 
puting a few additional trigonomet ric function. Usin g 
the direct angle relation, as in e.g. iBasak et al.l (j2QQ9h . 
may expose us to problems in the case of directions near 
the poles. The displacement field d(n) is obtained by 
taking the spherical gradient of a scalar potential, corre- 
sponding to the projected gravity field in the Born ap- 
proximation. 



Normalized error 

Figure 5. Accuracy of the error estimate - The thick black line 
sho ws t he measured error normalized using the error estimate, 
Eq. (|16|) , at each pixel. The overplotted thin red line is a Gaussian 
distribution of width 1 and centered on showing perfect agree- 
ment between predicted and actual error. 

3.2. Test of producing lensed CMB maps 

We now test our procedure for generating precise 
lensed CMB maps. We compare our generated maps to: 

- the true lensed maps, obtained by summing exactly 
the spherical harmonics at the position of interpo- 
lation 

- the naive interpolation procedure using a simple 
bilinear interpolation of CMB maps simulated at 
high resolution. This algorithm is included in the 
HealPix package for visualization purposes and 
was not intended for scientific use, but it is still 
useful as a point of comparison to assess whether a 
more complicated interpolation procedure is war- 
ranted. 

- the ECP bicubic interpolation algorithm imple- 
mented in LensPix. This is the current main- 
stream algorithm for quickly computing lensed 




T - T rei (K/K) 



xlO 



Q - Qrei(K/K) 



xlCT 



Figure 6. Error distribution of the Equicylindrical projection interpolation and of FLINTS - We represent the measured error distribution 
of the temperature (T, left panel) and one plane of polarization (Q, right panel) map. We use iVgide = 4, 096 and ^max = 4, 096 and no 
multiplication factor for the ECP method. The error distribution is computed by taking the difference of the value predicted by the 
interpolator to the exact value computed using ECP method. In black (red) solid line, we represent the error distribution of our Flints 
(ECP interpolation) algorithm. 



EC R T error map 



FLINTS, T error map 




an u 2 

Figure 7. Error distribution of the Equi- cylindrical projection interpolation and of FLINTS on the sky - We represent a comparison 
of the error distribution on the sky of the interpolated temperature (T, top panels) and one polarization plane (Q, bottom panels) map. 
The errors in the ECP interpolation, are represented in the left panels. The errors of the FLINTS interpolation are represented in the 
right panels. We compare to the exact lensing method. All maps were computed at Af side = 1,024, with Cax = 2,048. For visualization 
purposes, we degraded the error maps to a resolution of iV side = 128. 
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maps. 



As it is very expensive to compute the true lensed map 
on the full sky we limit ourselves to testing our method 
on a restricted subset of pixels distributed uniformly over 
the sky. More specifically we chose directions on the sky 
at a resolution of N s ^ e = 64 (49152 pixels). 

We show the result in Fig. 2J We note that Flints 
behaves much better than the naive interpolation. Our 
interpolation procedure, executed at a resolution of 
is able to match fairly well with the naiv e interpolation 
used at 4Af S id e - This matches iLewisI (j2QQ5h who indicated 
that one needs at least 16 times more pixels than the base 
CMB map to produce acceptable spectrum using a naive 
interpolation procedure. 

Moreover, the tails of the error distributions of the 
interpolated field are much more Gaussian, as illustrated 
in Fig. [5] There, we represented the error normalized by 
the expected standard deviation given by Eq. ([T6]) . If 
the errors are Gaussian with exactly this deviation, we 
must obtain a Gaussian of standard deviation equal to 
one, which is exactly what we obtain. 

In Table [2j we compare the performances of our 
method, F l ints, with the method described in 
iBasak et all (|2009h . labelled TORUS, and the ECP in- 
terpolation implemented in LensPix. For Flints and 
ECP we give an estimate of the attained precision. We 
also measure the time required to produce one lensed 
map. We note that the precision is better for Flints 
than for ECP. The problem for ECP are the heavy tails 
in the error distribution, as shown in Fig. [6] Further- 
more, at low A^ide fixed ^ max , the field fluctuations are 
undersampled which degrades the predictive properties 
of ECP method. On the other hand, Flints keeps er- 
rors lower because it takes into account the underlying 
fluctuations through the use of the angular power spec- 
tra in the weights. At N s ^ e = 4, 096, the central part of 
the error distribution represented in Fig. [6] is the same 
showing convergence of the two methods. 

A comparison of the differences in the spatial distribu- 
tion of the interpolation errors of the ECP and FLINTS 
method are given in Fig. [71 We represent there the sky 
distribution of the errors in the lensed temperature and 
polarization maps. We compute the reference maps using 
the full resummation of the spherical harmonics at the 
displaced positions. The maps are computed at a resolu- 
tion A^ide = 1,024, 4n ax = 2,048, with an oversampling 
factor equal to one for the ECP method. We note that 
the error distribution of the ECP interpolation is linked 
to the projection of the ECP mesh onto an HealPix 
mesh. On the other hand, the FLINTS interpolation is 
essentially tracking the shape of the HealPix grid, as 
already seen in Fig. 03 As in Fig. [6j the overall ampli- 
tude of error of the ECP interpolation is larger than the 
one given by the FLINTS interpolation. 

The total time for both methods is similar but Flints 
has initialization time overhead for precomputing S -1 
(Section l2.4p . This overhead may still be further opti- 
mized by using more of the symmetries of HealPix pix- 
elization. Overall, Flints is more precise than ECP by 
order of magnitude for simulations resembling data from 
a high resolution CMB polarization mission. Both of 
these methods are much faster than the TORUS method. 



3.3. Lensed power spectra 

As an additional test of the precision of our method we 
compare the power spectra of our lensed temperature and 
polarizat ion maps with the theoretical predictions from 
CAMB (jLewis et al]l2000UChallinor fc Lewisl l2QQ5). We 
show the difference of the average spectrum of 350 lensed 
maps and the unlensed spectrum in Figure [51 

We tested how our method fared compared to naive 
temperature interpolation and on a limited number of 
exactly interpolated pixels. We now check the inter- 
polation of both the temperature and the polarization 
fields on different scales by considering the difference 
between the lensed spectra, Cj^ nsed and the unlensed 
spectra. To assess the precision of our procedure, we 
use the spectra compute d by CAMB (Le wis et al.ll200Ci 
Challinor & Lewis 2005) as a reference. The results, 
given as the difference of the lensed spectra to the un- 
lensed spectra, are given in Figure [H We used A^de = 
4,096 and £ max = 5,000 as an input to FLINTS. The 
CAMB spectra were predicted using £ max = 10,000 and 
krj, max = 100,000 in very high accuracy mode. We also 
show the relative difference between our spectra and 
CAMB spectra. We see that the difference between 
the two does not exceed 0.5% statistically at t = 3,000. 
While the B-mode spectrum has a small but systematic 
excess of power of about 0.1 — 0.2% in our measured BB 
compared to the prediction given by CAMB, we note 
that this is of the same order as the advertised accuracy 
CAM B even in high-precision mode (jChallinor fc Lewisl 
l2QQ5h . In addition, CAMB is not guaranteed to be ac- 
curate at this level beyond i ~ 2000, so the comparison 
breaks down at this point and it is not clear whether 
CAMB or Flints is more accurate. 

In any case, the interpolation accuracy is more than 
sufficient for practical purposes in all cases. To illustrate, 
we display the cosmic variance range for the temperature 
power spectrum in the first panel of Figure [H This shows 
the unavoidable error in the estimation of Ci from a per- 
fect, all-sky, unlensed map. For the very small B-mode 
signal these numerical errors will remain much smaller 
than measurement error for the foreseeable future. 

4. CONCLUSION 

For a given interpolation stencil, we describe the op- 
timal interpolation technique for isotropic band-limited 
fields of arbitrary spin, sampled on the sphere. Taking 
advantage of the symmetry properties of the HealPix 
pixelization, the method is fast and memory-efficient. 
To test this approach we implement a supersampling fil- 
ter for HealPix temperature and polarization maps. A 
Monte Carlo study confirms both the predicted precision 
and our estimates of memory and CPU time scaling. 

We demonstrate this interpolation method to be pow- 
erful tool to simulate lensed CMB temperature and 
polarization maps from unlensed maps. Our Monte 
Carlo comparison to exact reference maps computed 
by LensPix and to predicted lensed power spectra 
by CAMB demonstrate that we achieve an accuracy 
which exceeds the requirements of the Planck data 
(jThe Planck Collaboration! |2006[ ) while reducing the re- 
quired computational time by an order of magnitude 
compared to using, e.g., the Torus method. In addi- 
tion, the method allows very easy parallelization as the 
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Figure 8. Precision of lensed power spectra - We show the differences between the lensed spectra and the unlensed spectra for ~350 
realizations of CMB fluctuations and lensing potentials. In the left panels, the dashed, red line is the difference as computed by CAMB, 
and the solid, blue line the result obtained using our interpolation technique. Cosmic variance error bars are shown as a dotted line in the 
top-left most panel. In the right panels, we represented the relative difference between the CAMB prediction and FLINTS prediction. We 
represented the change in the Cg for the temperature (top row), the E polarization mode (middle row), the B polarization mode (bottom 
row). The fluctuations are here normalized by the CMB temperature and are unitless. The red lines in the right panels represent perfect 
agreement with CAMB. The lensed maps were computed using a CMB map at N s ^ e = 4,096, Cax = 5,000 and 9 neighbors. 
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procedure is strictly local in pixel space. We compared 
the performance and the precision of our method to the 
Equi-cylindrical projection interpolation method. The 
two methods have similar speed. Flints is more precise 
and has no catastrophic errors. 

We conclude that this method is a very promising tech- 
nique in terms of speed, precision and scalability for the 
simulation of high resolution maps of the lensed CMB 
temperature and polarization anisotropics. Flints en- 
ables us to produce lensed maps as cheaply as making 
a spherical harmonic transform, and makes us capable 
of producing thousands of simulations of the lensed sky 
within an acceptable computational time. This advance 
may allow us to run a full likelihood analysis of the lens- 
ing potential in observed CMB data with current com- 
puter technology, which is not possible with other known 
methods of computing lensed maps. 
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APPENDIX 

A. SPIN- WEIGHTED SPHERICAL HARMONIC 
Spin s functions s f transform under a locally planar rotation R(6) about the direction h as 

J(R(6)n) = e is6 a f{h) (Al) 

On the sphere, these functions may be expanded on the spin- s spherical harmonic basis s Ye m (n) as (jNewman fc Penrose] 
[19661 : IGoldberg et IH[l967[ iZaldarriaga fc Selialdll997h 

+oo +e 

sf(n) = ^2 ^2 sfe y m[sYe,m(n)}. 

£=0 m=-£ 

As for spin-0 function, the spherical harmonic coefficient s fi,m ma y be obtained using 



(A2) 



dfi(n) [ s yj, m (n)]* s f(n) 



(A3) 
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with S 2 being the sphere in dimensio n three. The spin- weig hted spherical harmonic function may be expressed directly 
from the Wigner rotation matrices: (|Goldberg et 



spm-weigr. 
aD[l967|) 



V+ljf ( 



(A4) 



where we used the Condon-Shortley phase ([Condon fc Shortlevlll951l ) convention, and (j) are respectively colatitude 
and longitude on the sphere. The Wigner D matrix D £ m m , may be further expanded with the help of the Wigner d 
function 



D 



?,p) =e 



which yields 



2 ^ + 1 c im 
47T 



Furthermore, we recall the spin-s spherical harmonic addition relation (e.g. iHu fc White1ll997| ) 



[-'^m(» , ^ , )][-^m(fl,0)] = 



2£+ 1 

47T 
47T 

2£+l 
( 

47T 



&-s>A*'J^)&-sAhW) 



J, a) 



(A5) 
(A6) 

(A7) 
(A8) 

(A9) 
(A10) 



with (a, /?, 7) the Euler angles of the rotation bringing the direction (6 f ,<j) f ) to (#,0). 

B. COMPUTING THE WIGNER D FUNCTION 

We use t he decomposition of the Wigner d function in terms of their Fourier representation. This decomposition is 
taken from lEdmondsl (|1957| ). We start by factorizing a nodal rotation: 



R(p, P,0) = R (- J, 0, 0) R (0, ~,0)R 0, 0) R (o, |, o) R (|, 0, o) . 
Expressing this matrix multiplication in terms of the elements of the D matrices yields the identity: 



(Bl) 



(B2) 



If we let Af 



-mi,m 2 = d mi,m 2 (f ) and 5 n,m 1 ,m 2 = * mi m2 ^ mi Ai )m2j then we see explicitly the expression of d £ mum2 (f3) 
in terms of a discrete Fourier transform: 



u mi ,m 2 / v D n,mi,m2 



in (3 



(B3) 



The recursive formula for d-matrix can be adapted for the specific case of /? = 7r/2. This calculation yields the 
following recursion formula: 



with an initial condition 



A £ - 
^£,0 — 



21-1 
2i 



1/2 



x ^-l,0 



(£/2)(2£-l) 



1/2 



a£-1 

^-l,m 2 -l 



(£ + m 2 )(^ + m 2 - 1) 

^2 w 



{i-m 1 - l)(£ + m 1 + 2) 



(<-mi)(f + mi + 1 



1/2 



A^ 

mi+2,m 2 



A - 1 
^0,0 — l - 



(B4) 
(B5) 

(B6) 
(B7) 
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We start by using Eq. (|B4|) to increase the order in £ from the initial condition. At the same time, we compute all 
A £ £ m for < m < i using Eq. (|B5|) recursively. At the end of the recursion we have access to all A^ m at the required 
t. We use the last equation (|B6|) to compute A^ mi and A^ m2 for all — ^ < u < i at m\ and m2 fixed. 

The function m2 (P) is used to compute the angular correlation function described in Eq. (fl9|) . We decided to 
tabulate and interpolate using cubics the correlation function. For performance reason, we compute the Wigner-<i 
function at all f3 at sufficiently high resolution and then sum the contribution at each f3 for any given £ of the whole 
summation of d £ _ Q _ Q . 



